# define objects to save model estimates for plotting
diffs <- NA
diffs.se <- NA
beta1 <- beta2 <- NA
se1 <- se2 <- NA

VDEM <- read.csv("V-Dem-DS-CY-v6.1.csv")
VDEM <- subset(VDEM, select=c(COWcode, year, v2cltort, v2cltort_sd, v2clkill, v2clkill_sd, v2csreprss, v2csreprss_sd, country_name))
names(VDEM) <- c("COW", "YEAR", "v2cltort", "v2cltort_sd", "v2clkill", "v2clkill_sd", "v2csreprss", "v2csreprss_sd", "country_name")

year <- 1949:2010



# loop through all the different model specifications for this treaty variable
for(j in 1:length(year)){
    
    data <- VDEM
    
    data$VAR <- data$v2cltort
    data$VAR.SD <-data$v2cltort_sd
    #data$VAR <- data$v2clkill
    #data$VAR.SD <-data$v2clkill_sd
    #data$VAR <- data$v2csreprss
    #data$VAR.SD <-data$v2csreprss_sd

    OBS <- rev(cumsum(rev(as.numeric(table(data$YEAR)))))
    
    # load binary treaty variables
    treaty <- read.csv("hr_treaties.csv")
    treaty <- subset(treaty, select=c(ccode, year, cat_rn, cescr_rn, cerd_rn, ccpr_rn, cedaw_rn, crc_rn))
    names(treaty) <- c("COW", "YEAR", "cat", "cescr", "cerd", "ccpr", "cedaw", "crc")
    
    # this is not correct but necessary for the replications
    WRONG <- TRUE
    if(WRONG==TRUE)treaty$cat[is.na(treaty$cat)] <- 0
    if(WRONG==TRUE)treaty$cescr[is.na(treaty$cescr)] <- 0
    if(WRONG==TRUE)treaty$cerd[is.na(treaty$cerd)] <- 0
    if(WRONG==TRUE)treaty$ccpr[is.na(treaty$ccpr)] <- 0
    if(WRONG==TRUE)treaty$cedaw[is.na(treaty$cedaw)] <- 0
    if(WRONG==TRUE)treaty$crc[is.na(treaty$crc)] <- 0
    
    # load polity data
    polity <- read.csv("p4v2010.csv", na.strings = c("-66", "-77", "-88", "-99"))
    polity2 <- subset(polity, select=c(ccode, year, polity2))
    
    # load updated Gleditsch trade gdp and population data
    new.gdp <- read.delim("gdpv6.txt")
    new.gdp <- subset(new.gdp, select=c(statenum, year, rgdppc, pop))
    names(new.gdp) <- c("ccode", "year", "rgdpch", "pop")
    
    # merge the data together and make lags
    data <- merge(data, treaty, by.x=c("COW", "YEAR"), by.y=c("COW", "YEAR"), all.x=TRUE, all.y=FALSE, suffixes = c("",".treaty"))
    nrow(data)
    
    data <- merge(data, polity2, by.x=c("COW", "YEAR"), by.y=c("ccode", "year"), all.x=TRUE, all.y=FALSE)
    nrow(data)
    
    data <- merge(data, new.gdp, by.x=c("COW", "YEAR"), by.y=c("ccode", "year"), all.x=TRUE, all.y=FALSE)
    nrow(data)
    
    
    data.lag <- data
    data.lag$YEAR <- data.lag$YEAR + 1
    data <- merge(data, data.lag, by.x=c("COW", "YEAR"), by.y=c("COW", "YEAR"), all.x=TRUE, all.y=FALSE, suffixes = c("",".lag"))
    
    data <- subset(data, !is.na(v2cltort.lag) & !is.na(ccpr.lag) & !is.na(polity2.lag) & !is.na(rgdpch.lag) & !is.na(pop.lag) & YEAR>=year[j] & YEAR<=2010)
    nrow(data)
    
    newdata <- list()
    
    # take n draws from the posterior distribution and make n new datasets in the list object just created
    for(i in 1:50){
        newdata[[i]] <- data
        newdata[[i]]$draw.lag <- rnorm(cbind(rep(1,nrow(data))), mean=data$VAR.lag, sd=data$VAR.SD.lag)
        newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(data))), mean=data$VAR, sd=data$VAR.SD)
    }
    
    
    # define model formula for each specification
    #FML <- "VAR ~ draw.lag + ccpr.lag"
    #FML <- "VAR ~ draw.lag + ccpr.lag + polity2.lag"
    #FML <- "VAR ~ draw.lag + ccpr.lag + polity2.lag + log(rgdpch.lag)"
    FML <- "VAR ~ draw.lag + ccpr.lag + polity2.lag + log(rgdpch.lag) + log(pop.lag)"
    #FML <- "VAR ~ draw.lag + ccpr.lag + log(rgdpch.lag) + log(pop.lag)"
    #FML <- "VAR ~ draw.lag + ccpr.lag + log(rgdpch.lag)"
    #FML <- "VAR ~ draw.lag + ccpr.lag + log(pop.lag)"
    #FML <- "VAR ~ draw.lag + ccpr.lag + polity2.lag + log(pop.lag)"
    
    # define function to combine multe model estimates (this incorporates uncertatiny from the latent variables, see the supplementary appendix section F for more details about this function)
    milm <- function(fml, midata){
        xx <- terms(as.formula(fml))
        lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
        ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
        vcovs <- list()
        for(i in 1:length(midata)){
            tmp <- lm(formula=as.formula(fml), data=midata[[i]])
            lms[,i] <- tmp$coefficients
            ses[,i] <- sqrt(diag(vcov(tmp)))
            vcovs[[i]] <- vcov(tmp)
        }
        par.est <- apply(lms, 1, mean)
        se.within <- apply(ses, 1, mean)
        se.between <- apply(lms, 1, var)
        se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
        list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)
    }
    
    # call the milm function with the list of datasets newdata
    model <- milm(fml=FML, midata=newdata)
    
    # create object for printing during the loop and then print the results to screen
    temp <- data.frame(model$terms, model$beta, model$SE, model$beta/model$SE)
    
    print(paste("Model Output"))
    print(temp)
    
    print(as.numeric(temp[temp$model.terms=="ccpr.lag",2]))
    print(as.numeric(temp[temp$model.terms=="ccpr.lag",3]))
    
    
    
    #par(mar=c(5,7,5,5))
    
    VDEM.est <- as.numeric(temp[temp$model.terms=="ccpr.lag",2])
    VDEM.se <- as.numeric(temp[temp$model.terms=="ccpr.lag",3])
    
    
    
    beta1[j] <- VDEM.est
    se1[j] <- VDEM.se
    
}






# barplot of the differences across all 8 model specifications
par(mar=c(4,4,4,6))
par(mfrow=c(1,1))
barplot(beta1, space=0, ylim=c(-.35,.35), col=grey(.9), yaxt="n", border="white", xlim =c(1,60))
abline(h=0, col=grey(.5))
for(i in 1:7){abline(v=(c(2,12,22,32,42,52,62)-.5)[i], lwd=.75, col=grey(.75), lty=2)}
for(i in 1:length(year)){
    lines(c(i-.5,i-.5), c(beta1[i]+1.96*se1[i], beta1[i]-1.96*se1[i]), lwd=1)
    lines(c(i-.5,i-.5), c(beta1[i]+se1[i], beta1[i]-se1[i]), lwd=2)
}
#box()
#axis(side=1, at=c(.5:61.5), labels=rep("", 62))
axis(2, at=c(-.20, -.15, -.10, -.05, 0, .05, .10, .15, .20), las=2)
mtext(side=4, "Model Coefficient \nChanging Standard DV", line=3.0, cex=1.0)
mtext(side=3, "Human Rights Latent Treaty Variable Coefficients", line=-1.5, cex=1.25)

axis(side=1, at=c(.5:61.5), labels=rep("", 62))
axis(1, at=c(2,12,22,32,42,52,62)-.5, label=c(1950, 1960, 1970, 1980, 1990, 2000, 2010), tick=F, cex.axis=1.5)
axis(2, at=c(-.35,-.30, -.25, -.20, -.15, -.10, -.05, 0, .05, .10, .15, .20, .25, .30, .35), las=2)




